source("/scratch/jc227089/evo-dispersal/evostoch/evostochFunctions.R")
setwd("/scratch/jc227089/evostochData/")

fname<-"varh"
flist<-list.files(pattern=fname)
flist<-flist[!grepl(pattern="concat", flist)]
flist<-flist[!grepl(pattern="Allee", flist)]

nreps<-10

concat<-c()
#poplist<-vector(mode="list", length=nreps*length(flist))
for (ii in 1:length(flist)) {
	load(flist[ii])
	print(length(out))
	for (jj in 1:length(out)){
		if (is.null(out[[jj]]$pop)) next
		#poplist[[(ii-1)*nreps+(jj)]]<-out[[jj]]$pop
		pars<-rep(unlist(out[[jj]]$parameters), each=2)
		pars<-matrix(pars, nrow=2)
		if (dim(out[[jj]]$out)[1]==5) out[[jj]]$out<-cbind(rep(NA, 2), t(out[[jj]]$out))
		pars<-cbind(matrix(rep(c(pop=ii, rep=jj), each=2), nrow=2), pars, out[[jj]]$out)
		concat<-rbind(concat, pars)
	}
}
colnames(concat)<-c("pop", "rep", names(unlist(out[[jj]]$parameters)), 
		colnames(out[[jj]]$out))
colnames(concat)[16]<-'xlim'
concat<-cbind(concat, evol=concat[,"h2H"]>0, abs.xlim=abs(concat[,"xlim"]))

save(concat, file="concat_evostoch_h2walker.RData")


c.df<-as.data.frame(concat)
rm(concat)

c.no<-subset(c.df, c.df$evol==0)
c.yes<-subset(c.df, c.df$evol==1)
rm(c.df)

var.summ<-c()
for (jj in 1:nlevels(as.factor(c.yes$h2D))){
	for (ii in 1:nlevels(as.factor(c.yes$h2H))){
		c.evol<-subset(c.yes, as.numeric(as.factor(c.yes$h2H))==ii & as.numeric(as.factor(c.yes$h2D))==jj)
		c.no.evol<-subset(c.no, c.no$pop %in% c.evol$pop)
		
		
		dem.st<-summary(aov(abs.xlim~1, data=c.no.evol))[[1]]$'Sum Sq' #pure demographic stochasticity
		evol.st<-summary(aov(abs.xlim~Error(pop), data=c.evol)) # partitioned into within and between pop (between=founder, within = evol+demographic stoch)
		evol.st.win<-evol.st$'Error: Within'[[1]]$'Sum Sq' #evol+dem st
		evol.st.bw<-evol.st$'Error: pop'[[1]]$'Sum Sq' #initial founder event
	
		#proportion attributable to each effect
		dem.st.prop<-dem.st/(evol.st.win+evol.st.bw)
		found.st.prop<-evol.st.bw/(evol.st.win+evol.st.bw)
		evol.st.prop<-(evol.st.win-dem.st)/(evol.st.win+evol.st.bw)
		var.summ<-rbind(var.summ, c(h2H=as.numeric(levels(as.factor(c.yes$h2H)))[ii], h2D=as.numeric(levels(as.factor(c.yes$h2D)))[jj], round(c(demog=dem.st.prop, founder=found.st.prop, evoln=evol.st.prop)*100, 2)))
	}
}


####################################################
### Figures ###
var.summ[,3:5]<-var.summ[,3:5]/100
library(lattice)

zmat<-matrix(var.summ[,5], nrow=4)
xy<-seq(0.1, 1, 0.3)
pdf('Figures/h2walker.pdf')
	par(mar=c(5,5,2,2), cex.lab=1.4)
	#wireframe(var.summ[,5]~var.summ[,1]+var.summ[,2])
	contour(x=xy, y=xy, z=zmat, xlab=quote(italic({h^2}[w])), ylab=quote(italic({h^2}[d])))
dev.off()